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ABSTRACT: The organization of movement in a changing image provides a valuable source of 
information for analyzing the environment in terms of objects, their motion in space, and their 
three-dimensional structure. This movement may be represented by a two dimensional velocity field 
that assigns a direction and magnitude of velocity to elements in the image. This paper presents a 
method for computing the velocity field, with three main components. First, initial measurements 
of motion in the image take place at the location of significant intensity changes, which give rise 
to zero-crossings in die output of the convolution of die image with a V 2 G operator. The initial 
motion measurements provide the component of velocity in die direction perpendicular to die 
local orientation of the zero-crossing contours. Second, diese initial measurements arc integrated 
along contours to compute the two-dimensional velocity field. Third, an additional constraint of 
smoothness of the velocity field, based on the physical constraint dial surfaces are generally smooth, 
allows the computation of a unique velocity field. The details of an aigoridim are presented, with 
results of the algorithm applied to artificial and natural image sequences. 
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1. INTRODUCTION 


The organization of movement in a changing two-dimensional image provides a valuable 
source of information for analyzing the emironment in terms of objects, their motion in space, and 
their three-dimensional structure. It is not surprizing, therefore, that the analysis of motion plays a 
central role in biological vision systems, and in recent years, has played an increasingly important 
role in computer vision systems as well. For biological systems, the analysis of movement is crucial 
for such basic tasks as detecting and tracking prey, responding quickly to an approaching predator, 
and guiding locomotion through a complex emironment. Perhaps the most remarkable use of 
visual motion is the recovery of three-dimensional structure from the changing two-dimensional 
projection that a moving surface casts onto the eye. The ability of die human visual system to 
perform this recovery was first demonstrated in the studies of Wallach and O’Connell (1953) and 
Johansson (1973, 1975). For animals without binocular vision, motion is a primary cue to the 
three-dimensional structure of the world around them. 

The extensive use of motion by biological systems, and in particular the human visual system, 
demonstrates the feasibility of carrying out certain information processing tasks and helps to 
establish specific goals for the analysis of movement in time-varying imagery. This analysis divides 
naturally into two parts. The first is the measurement of motion; for example, the assignment 
of direction and magnitude of velocity to elements in the image, on the basis of the changing 
intensity pattern. The second is the use of motion measurements; for example, to separate the 
scene into distinct objects, and infer their three-dimensional structure. 

This paper presents a computational study of the measurement of visual motion. It is a 
problem that was found to be surprizingly difficult, both in computer vision, and in modelling 
biological vision systems. The main theoretical problems posed by this measurement are introduced 
in Section 2. Section 3 presents a particular formulation of the input and output representations 
for tire underlying computations. A theoretical analysis of the computation is then presented in 
Section 4, followed by the discussion of a specific algorithm in Section 5. Examples of the results 
of the algorithm suggest that the proposed computation is feasible for computer vision systems, 
and also yields solutions that are consistent with human motion perception. 

2. WHY IS THE MEASUREMENT OF MOTION DIFFICULT? 

The changing image may be represented by a two-dimensional array of time-dependent light 
intensities, I(x,y,t). Motion in the image may be described by a two-dimensional vector field 
\(x,y,t) that specifies the direction and magnitude of velocity at points with coordinates (x,y) at 
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Figure 1. A velocity field. The black lines represent local velocity vectors, displayed at evenly spaced 
points over the original image. 


time t. The measurement of visual motion may then be formulated as the computation o(\[x,y,t) 
from l(x,y,t). Figure 1 illustrates a velocity field that was computed from a sequence of aerial 
photographs, using the algorithm that is developed in this paper. In the background, a single image 
is shown with reduced contrast. The superimposed black lines represent local velocity vectors that 
describe the movement of features in the image, as the plane flics over the scene. For some visual 
tasks, it may be sufficient to compute only certain properties of the velocity field; for example, to 
respond quickly <o a moving object, motion must be detected, but not necessarily measured. Other 
tasks, such as the recovery of three-dimensional structure, require a more complete and accurate 
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Figure 2. The aperture problem, (a) An operation that views the moving edge E through the local aperture 
A can compute only the component of motion c in the direction perpendicular to the edge, (b) and (c) 
The perpendicular components of velocity for a translating circle and square. 

computation of the velocity field, as illustrated in Figure 1 (Ullman 1979a, 1983a, b; Clocksin 
1980; Prazdny 1980; Longuct-Higgins 1981; Longuet-Higgins and Prazdny 1981). 


The measurement of motion poses significant theoretical problems for a computational study. 
First, local motion measurements, obtained directly from the changing image, in general only 
provide one component of local velocity. This is a consequence of the aperture problem , illustrated 
in Figure 2(a) (Wallach 1976; Fcnncma and Thompson 1979; Horn and Schunck 1981; Marr and 
Ullman 1981; Adclson and Movshon 1982). If the motion of the edge F. is to be measured by a 
local motion detector that examines only an area A that is small compared to the overall extent of 
the edge, the only motion that can be extracted is the component c in the direction perpendicular 
to the local orientation of the edge. A local detector cannot distinguish between motions in the 
directions b, c and d in Figure 2(a). In Figures 2(b) and 2(c), a circle and square undergo pure 
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Figure 3. Ambiguity of the velocity field, (a) The arrows represent two possible velocity fields that are 
consistent with the changing image, (b) The curve Ci rotates, translates and deforms over time to yield the 
curve C 2 . The velocity of the point p is ambiguous. 

translation in the directions given by the vectors at the center of the figures. The vectors along the 
contours represent the local perpendicular components of velocity that can be obtained directly 
from the changing image. To compute the true motion of the figure, a second stage of analysis is 
required, that combines these local measurements. 

This combination stage faces a deeper theoretical problem, however; the movement of 
elements in the image is not determined uniquely by the pattern of changing intensities. Thus, 
the true velocity field is not determined uniquely from the initial local motion measurements. 
Two factors contribute to this ambiguity of motion. The first is the loss of information due to 
the projection of the three-dimensional world onto a two-dimensional image; multiple surfaces, 
undergoing different motions in space, may project to the same two-dimensional image. The 
second factor is the loss of information due to the projection into a pattern of changing intensity. 
The image that a surface projects onto the eye may not be sufficient to determine its movement 
in space. As an extreme example, a matte white sphere, rotating about a central axis, cannot be 
determined as such, on the basis of its projected image. 

Figure 3 presents two simple examples that illustrate the ambiguity of the velocity field. In 
Figure 3(a), the solid and dotted lines represent the image of a moving circle, at different instants 
of time. In the first frame (solid line), the circle lies parallel to the image plane, while in the 
second frame (dotted line), the circle is slanted in depth. One velocity field that is consistent with 
the two frames is derived from pure rotation of the circle about the central vertical axis, as shown 
to the left in Figure 3(a). (T he arrows represent a sample of the local velocity vectors along the 
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circle.) There could also be a component of rotation in the plane of the circle, about its center, as 
shown to the right in Figure 3(a). In addition, this changing image might represent the projection 
of a different three-dimensional curve that is deforming over time, giving rise to yet another 
projected velocity field. 1'his ambiguity is not peculiar to symmetric figures such as circles; it is a 
fundamental problem that is always present. In Figure 3(b). the curve (\ rotates, translates and 
deforms over time, to yield the curve ()•>. The motion of points from C\ to (■■> is again ambiguous 
(consider, for example, different possible velocities for the point />). In general, there arc infinitely 
many two-dimensional velocity fields that arc consistent with the changing image. 

To compute motion uniquely, additional constraint is therefore required, in the form of 
basic assumptions about the physical world that generally hold true. The main focus of this 
paper is the derivation of a particular constraint, the smoothness constraint , which relies on the 
physical assumption that surfaces are generally smooth, compared with their distance from the 
viewer. A smooth surface usually generates a smoothly varying velocity field when it moves. Thus, 
intuitively, we seek a velocity field that is consistent with the changing image, and varies smoothly. 
Unfortunately, die re arc still infinitely many possibilities. A single solution may be obtained, 
however, by finding the velocity field that exhibits the least amount of variation. In Section 4.3, 
this constraint is formulated more precisely, in a way diat yields a velocity field solution that is 
madiematically unique and physically plausible. 

3. THE INPUT AND OUTPUT REPRESENTATIONS 

The measurement of motion has been formulated as the computation of an instantaneous 
two-dimensional velocity field; such a formulation requires that motion in the image be roughly 
continuous. There arc alternate representations of visual motion that are not so restricted. For 
example, motion can be described by an explicit correspondence over time, between elements in the 
image that represent the same physical feature under motion (Ullman 1979a). Motion measurement 
in this case requires locating idcndfiable elements in the changing image, and matching diem 
over time. The input for a correspondence scheme may consist of a set of discrete frames, with 
large spatial separations between corresponding elements. r lhe perception of motion by the human 
visual system also does not require that a stimulus move continuously across the visual field. 
With appropriate spatial and temporal presentation parameters, a stimulus presented sequentially 
can produce die impression of smooth, uninterrupted motion, as in motion pictures (Wertheimer 
1912). Why, then, have we chosen a formulation of die motion measurement problem that relics 
on roughly continuous motion? 

Recent psychophysical investigation has suggested that in the human system, motion may 
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Figure 4. Initial processing of an image, (a) The original image, containing 320 x 320 picture elements, 
(b) The convolution of the image with a operator, (c) The resulting zero-crossing contours. 

be analyzed by two different systems, termed short range and long range processes by Braddick 
(1974, 1980), It has been proposed that the short range process analyzes continuous motion, or 
motion that is presented discretely, but with spatial displacements at most 10' - 15' of visual arc, 
and temporal intervals up to 60 - 100 milliseconds. The long range process would then analyze 
motion over larger spatial and temporal intervals. The ability of the human visual system to 
infer motion from the discrete displacement of image elements, over considerable distances and 
temporal intervals, suggests an underlying correspondence computation. Under these conditions, 
there is no continuous motion of elements across the image to be measured directly. Short range 
motion is roughly continuous, however; we propose that the measurement of motion by the short 
range process may be appropriately formulated as the computation of an instantaneous velocity 
field. 

With regard to the input representation for the velocity field computation, in order to detect 
movement in a changing image, there must be a variation of intensity over space and time; 
the combination of the two variations can be used to measure the direction and magnitude of 
velocity. The explicit comparison of spatial and temporal derivatives of intensity forms the basis of 
a class of motion measurement schemes referred to as gradient schemes (Limb and Murphy 1975; 
Cafforio and Rocca 1976, Fukinuki el at. 1976; Fennema and Thompson 1979; Horn and Schunck 
1981; Marr and Ullman 1981). In principle, motion measurements may be obtained wherever 
intensity varies. Marr and Ullman (1981) proposed, however, that initial motion measurements 
in the human system are made only at the locations of significant intensity changes. To detect 
these intensity changes, Marr and Hildreth (1980) proposed that a powerful operator for the initial 

7 





Figure 5. Decomposition of velocity, u- 1 ^) and u T (s) are unit direction vectors, perpendicular and tangent 
to the curve, and v x (s) and v T (s) are the two velocity components. 

filtering of an image is the Laplacian of a Gaussian, V 2 G (approximated in shape by the difference 
of two Gaussian functions). The elements in the output of an image convolved with V 2 G, which 
correspond to the locations of intensity changes, are the zero-crossings (Marr and Poggio 1979). 
Figure 4 shows an image that has been convolved with a V 2 G operator, and the resulting 
zero-crossing contours. Marr and Ullman (1981) proposed that initial motion measurements take 
place at the locations of these zero-crossings, using a mechanism that combines spatial and temporal 
gradients of the filtered image. Further motivation for this input representation can be found in 
Hildreth (1984a). 

In two dimensions, the initial measurements face the aperture problem. For the case of 
contours, local motion measurements provide only the component of motion in the direction 
perpendicular to the orientation of the contour. The component of velocity along the contour 
remains undetected. More formally, the two-dimensional velocity field along a contour may be 
described by the vector function V(a), where a denotes arclength. V(a) can be decomposed into 
components tangent and perpendicular to the contour, as illustrated in Figure 5. u T (a) and u- L («) 
are unit vectors in the directions tangent and perpendicular to the curve, and v T (a) and v4-(a) 
denote the two components: 


V(a) = v~*"(a)u~^(a) + «4-(a)u-*-(a). 


0 ) 


'I'hc component ^(a), and direction vectors u T (a) and u-Lf.,), arc given directly by the initial 
measurements from the changing image. I’hc component t> T (a) is not, and must be recovered. 



/"""A to compute V(fl). Intuitively, the set of measurements given by v 1 (.•<) over an extended contour 

should provide considerable constraint on the motion of the contour. Additional constraint is still 
required, however, to determine this motion uniquely. 

4. ADDITIONAL CONSTRAINT I OR MOTION MEASUREMENT 

in this section, the computation of die velocity' Reid is discussed from a theoretical viewpoint. 
The section is organized by the type of additional constraint that is used for this computation. 
Three types of constraint arc considered: (1) velocity is constant over an area of the image: (2) 
the velocity field is consistent with rigid rotation and translation of objects in the image plane; 
and (3) the velocity field is smooth, and exhibits the least variation among the set of velocity fields 
consistent with the changing image. 

The discussion of the constant velocity assumption, in Section 4.1, is primarily a review of 
previous work on motion measurement, as this assumption underlies most methods that have been 
proposed, both for computer and biological vision systems. The class of rigid motions in the image 
plane has been addressed in correspondence schemes for motion measurement. In Section 4.2, 
we present a scheme for analyzing the instantaneous velocity field for this class of motions. The 
„ constraint on smoothness of the velocity field, discussed in Section 4.3, was motivated by the work 

of Horn and Schunck (1981) on the optical flow computation, and allows the analysis of more 
general classes of motion. 

4.1 The Constant Velocity Constraint 

Much of the previous work in motion analysis has assumed that velocity is constant, at least 
over small areas of the image. If the projection of the scene onto the image plane is orthographic, 
the constant velocity constraint is strictly valid only when objects undergo pure translation. 
Examples of methods that use this constraint include those based on the cross-correlation of 
intensity, used both in computer vision (Lccse el al. 1970; Lillestrand 1972; Smith and Phillips 
1972; Wolferts 1974) and in modelling biological vision systems (Rcichardt 1961; Lappin and Bell 
1976; Petersik, Hicks and Pantle 1978; Anstis 1980), and correspondence schemes (Potter 1975, 
1977; Mulch and Thompson 1982; Lawton 1983). In addition, most gradient schemes assume that 
velocity is constant over an area of the image (Limb and Murphy 1975; Fennema and Thompson 
1979; Thompson and Barnard 1981; Marr and Ullman 1981). For gradient schemes, the constraint 
on velocity imposed by a single measurement of can be illustrated graphically in velocity 
space . a space in which the x and y axes represent the x and y components of velocity, which we 
denote by V x and V„. This is shown in Figure 6. When mapped to velocity space, the velocity 
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Figure 6. Velocity constraints in velocity space. t> x (.s) is the perpendicular component of velocity, and 
u- 1 ^) is the unit direction vector, at a point p on the image curve. The velocity vector at p must project 
to the line /; examples are shown with dotted lines. 

vector at a point on the contour must terminate along the line / perpendicular to the vector 
u J -(s)uT( s ); examples are shown by the dotted arrows. For the case of uniform translation, the 
lines of constraint formed by multiple measurements of tA(,,) along a contour intersect at a single 
point in velocity space, corresponding to the endpoint of the true velocity vector for the contour. 

Some gradient schemes for motion measurement make explicit use of this intersection point 
(Fennema and Thomspon 1979; Thompson and Barnard 1981: Adelson and Movshon 1982). 
Marr and Ullman (1981) proposed a zero-crossing based scheme, in which each local motion 
measurement restricts the true direction of velocity of a patch to lie within a 180° range of 
directions to one side of the zero-crossing contour. A set of measurements taken at different 
orientations along the contour further restrict the allowable velocity directions, until a single 
direction is obtained, which is consistent with all the local measurements. 

In many natural dynamic scenes, particularly those arising from motion of die observer, the 
velocity field can be approximated locally by pure translation. Most of die above schemes involve 
local operations for computing displacements or velocities, and can be effective for analyzing diis 
type of motion. In more general situations, for example, when a nearby surface rotates in depth 
and deforms, die assumption of local translation is no longer sufficient. 
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4.2 Rigid Motion in (lie linage Plane 


Some motion measurement schemes allow objects to undergo rigid rotation and translation 
in the image plane, l or example, Davis ct a!. (1983) present an iterative gradient scheme that 
combines motion constraints along contours, starting from points of known velocity, and utilizing 
the perpendicular components of velocity. Several correspondence schemes also allow the analysis 
of this class of motions (Ullman 1979a: Aggarwal and Duda 1975; Chow and Aggarwal 1977; 
Barnard and Thompson 1980; Nagel 1982). 


In the remainder of this section, we analyze die instantaneous velocity field for this class of 
motions, and present a simple geometric construction for computing die v elocity field. Suppose that 
a rigid curve undergoes general motion in space. Its instantaneous motion may be described as the 
combination of; (1) a rotauon with angular velocity ta about a single axis in space, which wc denote 
by the unit vector n = [n x ,n y ,n z } 7 ( T denotes the transpose of a vector), and (2) a translation, 
which wc denote by die vector t = [t x ,t y ,t z \ T . Let die curve be given by C = y(s), 2 (a)), 

where a denotes arclength. The location of a point on the curve may be given by the position 
vector r = [x(.s), y(s), •s(.‘*)j T . If it is assumed diat the axis of rotation passes through the origin 
of the coordinate system, the velocity of a point given by the position vector r is equal to the 
cross product of r with the vector <an. If wc dicn let die optical axis lie along the 2 -axis (with 
die positive 2 -axis directed toward the viewer), and let the projection of the curve onto the image 
plane (the (x,y) plane) be orthographic, then the two-dimensional velocity field V(«) along the 
contour is given by: 


V(a) = M(r x«n + t) = —u>z{s) 



-y[ s Y 

x(s) _ 



( 2 ) 


M denotes die matrix that performs the orthographic projection. The first term in the resuldng 
expression describes the component of die velocity field due to rotation in dcpdi about an axis 
parallel to the image plane (the axis n = [n I ,nj / ,0] r ); the second term is the component due 
to motion in the image plane (rotation about the axis n = [0,0, n.] r ), and the diird term is the 
translation component. 


Consider the restricted case of rigid motion in die image plane; the velocity field now 
corresponds to die combination of a translation, and rotation about the axis n = [0,0, l] r . Thus, 
V(a) is given by: 


V(«) 


-y(») 

x{s) 


+ 


A j 


(3) 


V(.s) is simply a translation, rotation and scaling of the image curve (*(«),?/(#)). as illustrated in 
figure 7. In Figure 7(a), the curve C\ undergoes a small rotation and translation in die image 




plane to yield die curve C 2 . The arrows indicate a sampling of the local velocity vectors along 
die curve. In Figure 7(b), these velocity vectors have been translated to a common origin in 
velocity space. The curve in velocity space has the same shape as the image curve Cj; its size is 
proportional to angular velocity w, and it is rotated 90° with respect to C t (this relationship is also 
used in kinematics (Hartenberg and Dcnavit 1964)). The addidonal translation of the curve in the 
image yields the same translation of the curve in velocity space. In the case of pure translation, 
the image of the velocity field in velocity space degenerates to a single vector. For the case of 
discrete motion of a curve, a similar relationship exists between the shape of the image curve, and 
its displacement field , the set of vectors that describe the discrete displacement of points on the 
curve (Hildreth 1984a). 


A simple geometric construction can be used to compute the velocity field, for the class of 
rigid motions in the image plane. If the true direction of velocity is known at two points on the 
contour, the direction of velocity can be computed everywhere as follows: (1) construct the lines 
perpendicular to the direction of velocity at the two known points, (2) compute die intersection 
of diese two lines, (3) from every point p,- along the contour, construct the line to die intersection 






Figure 8. Constructing the velocity field. The direction of velocity is known at points pi and p 3 , and is 
computed for the point p 2 . using the construction scheme described in the text. 

point; the true direction of velocity is perpendicular to this line. In Figure 8, the direction at P 2 is 
derived, given known directions at pi and p 3 . This construction is simply locating the point about 
which the motion can be described as pure rotation. For pure translation, the two lines, from 
points of known direction of velocity, arc parallel, so the direction of motion everywhere is equal 
to the direction of motion of tire known points. Certainly, if both the direction and magnitude 
of velocity were known at two points on the contour, the global motion parameters, and hence 
direction and magnitude of velocity everywhere, could be computed. From the direction of velocity 
alone at two points on the contour, however, the direction of velocity can be obtained everywhere. 
If, in addition, the perpendicular components of velocity are known, then both direction and 
magnitude of velocity can be computed along the contour. 

There are at least two sources for points of known velocity direction in the image. First, 
identifiable features, such as points of high curvature, may be tracked in two dimensions (Lawton 
1983). Second, points at which tire perpendicular component of velocity is zero are constrained 
to move along the direction of the tangent to the curve. For the case of a smooth, closed curve 
derived from a single object moving rigidly in the image plane, there exists at least two points on 
the curve for which the perpendicular component of velocity is zero (Hildreth 1984a). Suppose 
the velocity field computation is focused at the locations of zero-crossings derived from the image. 
Zero-crossing contours arc generally closed (except at image boundaries), so there is usually 
sufficient constraint from the image to compute the velocity field, for this simple class of motions. 

The methods described in this section, for computing the velocity field for the restricted 
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class of rigid motion in the image plane, arc not sufficient for general motion analysis; however, 
they may be useful for the initial detection and rough measurement of motion in the periphery, 
the analysis of motion during smooth pursuit eye movements, or the recovery of observer motion 
from optical flow (Pra/dny 1980; I .onguet-I liggins and Pra/dny 1981). In computer vision, there 
arc restricted applications for these techniques, such as the tracking of objects, or computation of 
camera motion (Briiss and Horn 1983; Lawton 1983). 

4.3 l'he Smoothness Constraint 

In this section we derive a more general constraint on the velocity field, that allows the 
computation of the projected motion of three-dimensional surfaces that move freely in space, and 
deform over time. We rely on the physical assumption that the real world consists predominantly 
of solid objects, whose surfaces arc generally smooth compared with their distance from the viewer. 
A smooth surface in motion usually generates a smoothly varying velocity field. Thus, intuitively, 
we seek a velocity field that is consistent with the motion measurements derived from the changing 
image, and which varies smoothly. Unfortunately, there is an infinity of velocity fields that satisfy 
these two properties. Horn and Schunck (1981), in their work on the optical flow computation, 
suggest that a single solution may be obtained by finding the velocity field that varies as little as 
possible. In the remainder of this section, we show how this constraint may be formulated more 
precisely, in a way that guarantees a velocity field solution that is mathematically unique and 
physically plausible. 

4.3.1 Measuring Variation in Velocity 

To find the velocity field that varies the least, some means of measuring the variation in 
velocity along a contour is required. This can be accomplished in many ways. For example, we 
could measure the change in direction of velocity as we trace along the contour. Total variation in 
velocity could then be defined as the total change in direction over the entire contour. A second 
possibility is to measure the change in magnitude of velocity along the contour. Third, die change 
in the full velocity vector could be measured, incorporating both the direction and magnitude of 
velocity. Other measures are also possible. The goal of die computation is to find a velocity field 
that is consistent with die changing image, and minimizes one of these measures of variation in 
velocity. 

A measure of variation may be described more formally by defining a mathematical functional, 
0, which maps the space of all possible vector fields (along die contour), V, into the real numbers: 
0:V h-* i ft. This functional should be such that the smaller the variation in the velocity field, the 




c. d. 

Figure 9. Measuring variation in velocity, (a) The vectors V(«) are displayed at two nearby points on 
the image curve C. (b) The velocity vectors drawn in velocity space, where is indicated by the dotted 
arrow, (c) The direction of velocity for points on the contour is represented by the angle tp. (d) The velocity 
vectors of (c) are drawn in velocity space, where is shown. 

smaller the real number assigned to it. Two candidate velocity fields may then be compared, 
by comparing their corresponding real numbers. This formulation allows the development of an 
explicit method for computing the velocity field of least variation. A set of functionals can now 
be derived, based on the measures of variation that were previously mentioned informally: (1) 
variation in the full velocity vector, V(a), (2) variation in the direction of velocity, and (3) variation 
in the magnitude of velocity, all with respect to the curve. 
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1. \ ariat ion in V(,s) 

The local change in V(.s) with respect to the contour is given by the vector A scalar 
measure may be obtained by taking its magnitude: . In Figure 9(a). two nearby velocity 

vectors along the image contour C are shown. The vectors are translated to a common origin in 
velocity space in Figure 9(b). where the vector is shown as a dotted arrow. A measure of the 
total variation in the velocity field over an entire contour may be derived by integrating this local 
measure, suggesting a functional such as: 


0(V) 


/ 


ax 




da. 


Variations on this functional may also be considered, involving higher order derivatives, or 
higher powers, such as: 


0(V) 


/ 


d 2 X 

da 2 


ds 


or 


«V>- S\l 


ds. 


r^ 


2. Variation in Direction 

Let the direction of velocity be described by the angle <p, measured in the counterclockwise 
direction from the horizontal, as shown in Figure 9(c). In Figure 9(d), the local change in direction 
for two nearby velocity vectors along the image contour, given by , is shown in velocity space. 
Total variation in direction along the contour may again be obtained by integrating this local 
measure, leading to functionals such as the following: 



or variations involving higher order derivatives, or higher powers. 


3. Variation in Magnitude 

Finally, the total change in magnitude of velocity alone could be measured, using functionals 
such as: 


Again, we could also consider variations on this measure. 
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The functional that is used to measure variation may also incorporate a measure of the 
velocity field itself, rather than strictly utilizing changes in the velocity field along the contour, 
by incorporating a term that is a function of |Vj. This might he useful if we sought a velocity 
field that also exhibits the least total motion. In addition, the functional could become arbitrarily 
complex in its combination of |^|, or higher order derivatives. 

Given that there are many possible measures of variation, what criteria can be used to 
choose a single measure? First, from a mathematical point of view, there should exist a unique 
velocity field that minimizes the particular measure of variation; this requirement imposes a set 
of mathematical constraints on die functional. Second, the velocity field computation should yield 
solutions that arc physically plausible. These two criteria are important for die evaluation of 
any assumption diat is proposed for the motion measurement computation. If, in addition, we 
suggest that such a constraint underlies die motion computation in die human visual system, the 
minimization of this measure of variation should yield a velocity field that is consistent with human 
motion perception. 

4.3.2 Mathematical Uniqueness of the Velocity Field 

An examination of these measures of variation from a mathematical viewpoint suggests that 
a measure incorporating the change in the full velocity vector is necessary for the velocity field 
computation. The use of functionals that incorporate only a measure of direction or magnitude of 
velocity, for example, does not in general lead to a unique velocity field solution. It can be shown, 
however, that given a simple condition on the constraints that are derived from the image, there 
exists a unique velocity field that satisfies these constraints, and minimizes die particular measure 
of variation given by: / |^| 2 ds. The basic mathematical quesdon is, what conditions on the form 
of the functional, and die structure of the space of velocity fields, are needed to guarantee the 
existence of a unique solution? These conditions are captured by die following dicorem from 
functional analysis (Rudin 1973); 

Theorem : Suppose there exists a complete semi-norm 0 on a space of functions II, and 
that 0 satisfies the parallelogram law. Then, every nonempty closed convex set E C II 
contains a unique element v of minimal norm, up to an element of the null space. Thus, 
the family of minimal functions is 

{v +«| s e 5} 

where 

S — {w | v + w £ E} n U 


n 




and M is the null space of the functional 



M — {« | (-)(«) —- ()}. 

The functional (-) - {/ 1is a complete semi-norm that satisfies the parallelogram law, and 

the space of possible velocity fields that satisfy the constraints derived from the changing image, 
is convex (Hildreth 1984a, b) It then follows from the above theorem that this space contains a 
unique element of minimal norm, up to possibly an element of the null space. The null space in 
this ease is the set of constant velocity fields, because the condition that /] -T/.s = 0 implies that 

^7 = 0 everywhere, which implies that V(») is constant. In addition, the smoothness measure 
is non-negative, so that minimizing {/ is equivalent to minimizing /|^| 2 <fa. These 

further observations lead to die following result (Hildreth 1984a, b): 

Ifv-i-(s) is known everywhere along the contour, and there exists at least two points at 

which the local orientation of the contour is different, then there exists a unique velocity 
field that satisfies the known velocity constraints and minimizes / \ i3 j~\ 2 ds. 

An extended straight line docs not yield measurements at two different orientations, but in all other 
eases, there is sufficient information along a contour to guarantee a unique velocity field solution. 
The smoothness constraint can be used to compute a projected two-dimensional velocity field for 
any three-dimensional surface, whether rigid or nonrigid, undergoing general motion in space. 
While it is not yet clear whether this general formulation of the smoothness constraint, or die 
particular measure J\^-\ 2 ds, is the most appropriate for the motion computation, it is important 
that this measure satisfies certain essential mathemadcal requirements, that the other measures do 
not. The use of a functional incorporating only a measure of velocity direction, for example, which 
attempts to make die local velocity vectors as parallel as possible, docs not lead to a unique velocity 
field solution. In Appendix A.l, it is shown that for the simple case of a line segment moving 
rigidly in space, with known velocity vectors at its endpoints, diere does not exist a unique velocity 
field that minimizes either variation in direction, or variation in magnitude of velocity. This result 
is not surprising; velocity has two degrees of freedom, and can be decomposed in various ways (for 
example, into direction and magnitude, x and y components, normal and tangential components, 
or radial and rotational components). If the measure of variation in velocity constrains only one 
of these degrees of freedom, such as direction, the other is allowed to vary freely, and cannot be 
specified uniquely. The simpler functional, /|^-|d.i, also docs not yield a unique velocity field in 
general. In Appendix A.2, an example is presented for which a velocity field that minimizes this 
measure docs not exist. Measures involving higher order derivatives impose a higher degree of 
smoothness on the velocity field, which may not be necessary for this computation. Our approach 
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is to consider first the simplest measure, /|'}' |-< 1». If this is not sufficient for the velocity field 
computation, we can then consider a more complex measure. 


/■"N 


4.3.3 Physical Plausibility of the Velocity Held Solution 

The second criterion for evaluating a particular measure of variation in velocity is the physical 
plausibility of the resulting solution. One question that can be asked is. under what conditions 
will the velocity field that minimizes JVg-l 2 ^ be the correct physical velocity field? If we assume 
orthographic projection of the scene onto the image, there are at least two classes of motion for 
which this is true. The first consists of arbitrary rigid objects undergoing pure translation. In this 


case, fi 


0 everywhere along contours in the image, and hence jT/.i = 0. Since zero is 
the smallest value that tire measure can obtain, it follows that if there exists a valid solution that 
is consistent with pure translation, then this solution minimizes /ifjl 2 ^. Consequently, motion 
measurement schemes that rely on pure translation address a special case of this more general 
method. 

The second class of motions includes rigid polyhedra, undergoing general motion in space. 
In particular, the following can be shown (Hildreth 1984a, b): 

Suppose that a rigid three-dimensional object, consisting of straight lines intersecting 
in space, projects onto the image plane, using orthographic projection, in such a way 
that line intersections are preserved (that is, two lines intersect in the image if and 
only if their generators intersect in space). Further, suppose that this object undergoes a 
general displacement in space. Then the two-dimensional velocity field that satisfies tC-(a) 
measured along lines in the image, and minimizes J\‘)X\ 2 ds, is the correct projected 
two-dimensional velocity field. 

The proof of this result takes advantage of the fact that die velocity field varies linearly along 
a straight line that is moving rigidly in space. The classes of blocks world objects (Roach and 
Aggarwal 1979), and simple polygons in the image plane (Aggarwal and Duda 1975), are special 
cases of diis class of motions, for which an algorithm that minimizes /j^-| 2 d.s is guaranteed to 
compute the correct velocity field. 

Recently, Yuillc (1983) derived a general condition under which the velocity field that 
minimizes /|^| 2 ds is die correct velocity field. Let V'(.i) denote the true projected two-dimensional 
velocity field for a curve in motion, and let T(.<) denote the tangent vector along die curve. If die 
following relationship holds at every point on the curve: 
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then the velocity field V(.s) that satisfies the constraints imposed by r ! (,s) and minimizes /|^ |'-b/.s 
is the true velocity field V'(.s). The two classes of motion mentioned above correspond to eases for 
which — 0 along the curve, so that this general condition holds trivially. Yuille (1983) notes 
that this condition is also satisfied when a curve undergoes pure expansion in the image plane. 

For the class of smooth curves, undergoing general motion in space, we rely on an empirical 
investigation for demonstrating the physical plausibility of die velocity field solution. In general, 
the velocity field of least variation, in this case, is not the physically correct one. It is often 
qualitatively similar to the true velocity field, however. As we will show in Section 5.4,1, when 
die two velocity fields differ significantly, it appears diat die smoothest velocity field may be more 
consistent with human motion perception. 

A possible constraint on die velocity field, that is not considered explicitly, is the rigidity 
of the underlying surface. The computation of a smoothest velocity field does not necessarily 
seek a solution that corresponds to rigid motion, in either two or diree dimensions. This may 
at first seem physically implausible. When a three-dimensional curve rotates in space, however, 
its two-dimensional projection may undergo significant distortion in die image. Without knowing 
the three-dimensional structure of the curve, it is very difficult, if not impossible, to find a 
two-dimensional velocity field that corresponds to a single rigid motion in three dimensions. 
It is also the case that some of the motion that we encounter arises from surfaces that are 
nonrigid. If the analysis of motion is a two-stage process, with the measurement of two-dimensional 
motion preceding the derivation of three-dimensional structure from motion, a constraint such as 
smoothness may be the most restrictive type of constraint that may be used, which yields a unique 
solution, and still allows the analysis of general motion. 

5. AN ALGORITHM AND EXAMPLES 

In the previous section, die use of the smoothness constraint led to the formulation of die 
velocity field computation as an optimization problem. We seek a velocity field solution diat 
satisfies the constraints derived from the changing image, and minimizes the measure of variation 
along contours given by: f\^-\ 2 da. The computation may also be described as seeking a solution 
for which neighboring velocity vectors are as similar as possible. To further test the adequacy of 
this approach, it is necessary to specify an algorithm diat embodies the smoothness constraint, 
and examine the results of die algorithm for a number of motion sequences. An algorithmic 
study first allows us to demonstrate the biological feasibility of this approach, by showing that 









this formulation of the velocity field computation lends itself naturally to algorithms that involve 
simple, local, parallel operations (Ullman 1979b; Grimson 1981; Marr 1982). Second, an empirical 
study provides a further means for testing the physical plausibility of the resulting velocity field 
solution. 

In general, there are many algorithms that implement a given computational theory. For the 
velocity field computation, there arc many choices for the design of an algorithm that embodies the 
additional smoothness constraint. In this section and in Appendix B, we use a standard, iterative 
algorithm from mathematical programming, known as die conjugate gradient algorithm (I tienbergcr 
1973), to implement this computation. This particular algorithm is certainly not appropriate as a 
model for human vision. Our aim is to test the basic idea of computing the velocity field of least 
variation. If the results of the algorithm support the feasibility of this idea, we can then explore 
alternative algorithms to implement the theory, that may be more appropriate for the human 
system. 


5.1 The Discrete Formulation 


In Section 3.3, a continuous formulation of the smoothness constraint was presented; the 
image, however, is discrete, and therefore image contours consist of a set of discrete points. The 
first step in the design of an algorithm is to convert the continuous formulation into a discrete 
one. To do this, the functional can first be expressed in terms of the x and y components of* 
velocity, which we denote by \ z and \ y . The continuous functional becomes: 



(5) 


Suppose there are n points along a contour, whose coordinates are given by: 


{(zi> Vi ), { x 2, 2 / 2 ), (*n, 2/n)}• 


Let (V x ., \ y .) denote the velocity at the point (aAssuming for now that the contour is closed, 
and the points are evenly spaced along the contour, the following function defines a discrete 
correlate to the continuous functional: 

*,-± [('V - v,,_, r + (V„ - V,,.,)* ] + [(V„ - V,.)* + <V„ - V,.)*]. (6) 

i=2 

For the case of an open contour, tire discrete differences at the endpoints change. If points on 
the contour arc not evenly spaced, the discrete differences must be divided by the distances 
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between adjacent points on the contour. The goal of the computation is to compute the set of x 
and y components of velocity, {(V*,, V„,), (V*.,, \ v . { ),..., (V Xn , V„.)}, which minimizes the discrete 
function <t>|. 

5.2 Satisfying the linage C onstraints 

The second step toward' defining an algorithm is to incorporate the constraints on velocity 
that are derived from the changing image. Here there is a choice; the velocity field can be forced 
to satisfy the constraints exactly, or satisfy them approximately. In the first ease, the constraints 
can be expressed as follows: 

V ■ u 1 - D- 1 - = 0. (7) 


The equation states explicitly that the perpendicular component of velocity for the computed 
velocity field, V • u- 1 -, is equal to the measured perpendicular component, tk. (For simplicity, 
we omit the argument, a, to these functions, throughout the remainder of the paper.) Letting 
V = (V XI V,) and u- 1 - = (uf ,uf ), the constraints are simple linear constraints: 

V* v uf + V p< UjL - of = 0. (8) 


This equation can be used to express the y components of velocity in terms of the x components: 


' v. 


Ujs 


(9) 


This expression is valid for uf ^0. If uf = o, then uf = 1, and it follows that V Zi = „f , and 
V„. is unknown. This situation arises at points on the contour for which the vector u- 1 - is oriented 
horizontally (the contour is oriented vertically). At such a point, tk specifies the x component of 
velocity directly, and the y component is unconstrained. For points at which uf ^0, the above 
expressions for \ yi can be substituted directly into the discrete function At points for which 
uf = 0, the known values of V Xj and unknown parameters V y< can be substituted into <[>!. 
Given n points on the contour, there remains n unknown parameters to compute, with no further 
constraint on these parameters. 


In general, there will be error in the measurements of o-k From a practical standpoint, 
it may be advantageous to require that the velocity field only approximately satisfy these image 
constraints. This can be accomplished by requiring that the difference between V • uT and the 
measured v- 1 be small, rather titan requiring this difference to be zero, as in Equation (7). The 

continuous functional can be extended as follows: 

, 2 / \ 2l 


e “ / (^f) + (lv) *+ />/['' 


( 10 ) 
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p is a weighting factor that expresses the confidence in the measured velocity constraints. In this 
particular formulation, ft is constant over the entire contour, but in principle, p could vary, if 
confidence in the precision of individual measurements of <> 1 varies along the contour. The second 
term describes the least squares difference between the computed and measured perpendicular 
components of velocity. There will also be error in the measurement of the direction vector u L ; 
some of this error is captured implicitly in the second term of the above functional. The continuous 
functional leads to the following discrete function: 


c|» a 




pYI + 'W; 


(ii) 


The goal of die computation in this case is to find the set of x and y components of velocity, 
{(V Xl , \ 7 J/t ), , V„ 2 ),(V*„, Vj,J}, which minimizes the discrete function 4> 2 . These components 

satisfy the system of 2 n linear equations given by: 


< 94>2 

dV7 


= o, 


d^2 


— 0 , 


1 < i < n. 


( 12 ) 


In general, there may be several hundred points on a given contour, giving rise to a system 
of several hundred linear equations. To solve this system of equations, we use techniques from 
mathematical programming. 


5.3 An Algorithm from Mathematical Programming 


The general mathematical programming problem can be stated as: 
minimize /(x) 

subject to hi{\) = 0 i' = l,..., m 
ffy(x)<0 j = 1, .... r 
x e S 

where x is an n-dimensional vector of unknowns, x = (x 1( x 2 , ...,x n ). The objective function /, 
and constraints h t , i — and gj, j — l,...,r are real-valued functions of the variables 

xi,x 2 ,...,x n . The set S is a subset of the n-dimensional space. A particular optimization problem 
may or may not give rise to constraints of tire form A,(x) = 0 or gj(x)<0; that is, the problem 
may be a constrained or unconstrained optimization problem. 

For die case in which die image constraints are satisfied approximately, die velocity field 
computation can be formulated as an unconstrained optimization problem. In this case, we let: 


x = {V„, V, a , V >a> 
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The objective function /(x) -- 4>o. We seek a vector \ that minimizes /(x). with no further 
constraints on x. 

I'or the particular examples presented in this paper, the only constraints obtained from the 
image are the measurements of a 1 along contours. In general, there may be points at which V(.s), 
or the direction of velocity alone, is known. This information may be the result of a process that 
tracks localizablc features in the image, such as corners, isolated points, or line endpoints. These 
additional constraints may be incorporated easily into the velocity field algorithm (Hildreth 1984a). 

There arc several standard algorithms from mathematical programming that can be used to 
solve this particular optimization problem (I .ucnbcrger 1973). For the examples in the next section, 
die conjugate gradient algorithm was used to obtain the solution. This is an iterative algorithm 
that utilizes the gradient of the objective function /(x) to choose an optimal path to follow toward 
the final solution. The initial velocity field is given by the perpendicular velocity vectors, v-J-u- 1 -, 
along a contour. If there are n parameters to be computed, the conjugate gradient algorithm is 
guaranteed to converge to the exact solution in at most n iterations. Details of the application of 
the conjugate gradient algorithm to the velocity field computation are given in Appendix B. 

5.4 Examples of the Smoothest Velocity Field 

In this section, we present examples of die velocity field of least variation, for a variety - 
of motion sequences. The conjugate gradient algorithm was chosen for its fast convergence and 
because it is guaranteed to compute die velocity field of least variation. The algorithm provides 
us with a means of testing whether the use of the smoothness constraint is a reasonable approach 
to the velocity field computation. 

5.4.1 Ideal Smooth Curves 

We begin with some simple curves, undergoing rigid motion. For diis first set of examples, 
the curves and dieir perpendicular components of velocity were generated analytically from a 
known velocity field, and therefore represent ideal input for die algoridim. Many of the examples 
were chosen because perceptual studies indicate that human observers sec motions that differ from 
die true motion of die curves. 

/. Rigid translation in the image plane 

In Figure 10(a), a sampling of die true velocity field is shown for a simple polygon, translating 
rigidly in the image plane. Figure 10(b) illustrates the perpendicular velocity vectors v-i-ul, which 
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form the input to die algoridim. In diis case, the velocity field of least variation, shown in Figure 
10(c), is identical to the true velocity field, as expected. 

2. Polygon Rotating in the Image Plane 

Figure 11 illustrates the true velocity field, initial perpendicular velocity vectors, and smoothest 
velocity field for a polygon rotating rigidly in the image plane, about the point 0. From the 
theoretical results of Section 4.3.3, we expect to obtain the correct velocity field in this case. Figure 
11(c) shows that the computed smoothest velocity field is in fact the correct one. 

3. Rotating Ellipses 

Figure 12 illustrates the true velocity field, initial perpendicular velocity vectors, and smoothest 
velocity field for an ellipse rotating rigidly in the image plane, about the point 0. In this case, the 
smoothest velocity field is quite different from the true velocity field. There is a reduced rotational 
component of velocity, and added radial component in the computed velocity field. The difference 
in total variation is significant. 

At first glance, one might not consider the smoothest velocity field in this case to be a 
plausible solution. In some earlier perceptual experiments by Wallach ct al. (1956), however, they 
noted that a rigid ellipse docs not appear rigid under rotation; it appears to deform continuously. 
In their experiments, simple geometric figures were placed on a rotating turntable, and observers 
described the perceived motion of the figures while fixating die center of the turntable (conditions 




of free and tracking eye movements were also used). Ellipses of various aspect ratios were observed. 
It was found that when an ellipse whose axes measured 25 and 23.5 cm was rotated about its 
center, it appeared to stand still while its contour pulsated. The largest effects were observed 
for an ellipse whose aspect ratio was 3:2; the entire figure appeared fluid, undergoing a strong 
deformation, as well as a rotation. For some observers, the deformation was more restricted, and 
did not occur in the immediate vicinity of the poles of the major axis of the ellipse. As the aspect 
ratio of the ellipse was increased, it appeared more rigid. The perceived deformation of the ellipse 
was the same, regardless of whether it rotated about its center, or was placed eccentrically on the 
turntable. 

For the ellipse of Figure 12, the aspect ratio is 2:1. The computed velocity field of least 
variation clearly implies a significant distortion of the contour, in addition to a rotation. In the 
immediate vicinity of the poles of the major axes, the smoothest velocity field is very similar to 
the true velocity field. In Figure 13, the true and smoothest velocity fields are shown for rotating 
ellipses whose aspect ratios arc 25:23.5 and 5:1, When the ellipse is nearly circular, the smoothest 
velocity field indicates a strong inward and outward pulsation of the contour (Figure 13(b)). The 
smoothest velocity field for the narrower ellipse, showm in Figure 13(d), is closer to the true 
velocity field than in the case of the ellipse with aspect ratio 2:1, implying less distortion of the 
contour. Finally, the smoothest velocity field for an eccentrically rotated ellipse differs from that 









Figure 12. Rotating ellipse, (a) The true velocity field for an ellipse rotating rigidly in the image about 
/'"'y the point O. (b) The initial perpendicular velocity vectors, (c) The computed velocity field. 

of tite centrally rotated ellipse only by the addition of a uniform translation along the contour. 
Thus, the same deformations of the ellipses is implied by the smoothest velocity field obtained 
for eccentric rotations. We conclude that the perception of the movement of rotating ellipses is at 
least qualitatively consistent with the computation of the velocity field of least variation. 

4. Rotating Spiral 

It is well known that a spiral appears to expand or contract, when it undergoes pure rotation 
about its center (Holland 1965). The perceived velocity field thus contains a large radial component, 
while the true velocity field contains only a rotational component of velocity. Figure 14 illustrates 
the tme velocity field, initial perpendicular velocity vectors, and computed smoothest velocity 
field, for a single arm of a rotating logarithmic spiral. The smoothest velocity field exhibits a large 
radial component of motion at the center of the spiral, which decreases toward the periphery. 
The perception of the movement of the spiral is qualitatively more consistent with the smoothest 
velocity field, than with cither the true or initial velocity fields. 
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Figure 14. Rotating spiral, (a) The true velocity field for a logarithmic spiral rotating in the image about 
the point O. (b) The initial perpendicular velocity vectors, (c) The computed velocity field. 

horizontal. As viewing distance increased (displacements decreased), the apparent direction of 
motion became increasingly more vertical. 


The predictions of the velocity field algorithm are shown in Figures 15(b) through 15(e). 
In this case, the two images were convolved with a V 2 G operator, and the zero-crossings 
(which encircle individual dots) were obtained. The time derivative, §- t [V 2 G*l) was computed by 
subtracting the two filtered images. At the location of a zero-crossing, the local vector u- 1 was 
computed from the gradient of the filtered image, and the perpendicular components of velocity, 
v-h, were computed as follows: 


t;-L — 



(13) 


where V denotes the gradient operator, and I' is the filtered image. In the leftmost column, the 
actual displacements between dots arc given as a function of w, the diameter of the central positive 
region of the initial V 2 C operator. The resulting smoothest velocity fields are shown enlarged in 
the rightmost column of Figures 15(b) through 15(c). As the absolute displacement of the dots 
decreases, the predicted direction of motion becomes increasingly more vertical, consistent with 
the perception. 
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Figure 15. The moving dot demonstration, (a) Gosed circles represent a dot from the first frame, and 
open dots represent dots from the second. The double arrows show the perceived direction of motion of 
the full dot pattern. The inserted scale indicates the size of the displacements, (b) through (e) The leftmost 
column indicates displacements as a function of w, the middle column illustrates the initial perpendicular 
components of velocity, and the rightmost column illustrates the computed velocity fields. 
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Figure 16. The split herringbone demonstration, (a) Alternating columns of a herringbone pattern of lines 
move in opposite directions, (b) and (c) The zero-crossing contours derived for two different size V 2 G 
operators, (d) and (e) The initial perpendicular components of velocity, (f) and (g) The computed velocity 
fields. 
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7. The "split herringbone" demonstration 

Adelson and Movshon (1983) present a motion demonstration, referred to as the "split 
herringbone." The display, shown in Figure 16(a), consists of alternating columns of line segments 
that tilt left or right in the odd and even columns. The odd (left-tilting) columns move upward, 
while the even (right-lilting) columns move downward. There arc two perceptions of motion that 
result, depending on the viewing conditions. When the display is of high contrast, sharp, and 
centrally fixated, the columns arc seen mming vertically, in opposite directions. When the display 
is of low contrast, blurred with a diffusion screen, or viewed peripherally, an illusion of rightward 
motion is seen. To explain this phenomenon, Adelson and Movshon (1983) distinguish two kinds 
of velocity field computations. The first is a local computation, that involves die tracking of 
localizablc features, such as the line endpoints in the herringbone pattern. The second is a more 
"global" process that combines velocity constraints derived from different parts of the image. They 
argue that the perception of vertically moving columns represents the results of the first type of 
process, while the perception of rightward movement represents the results of the second. 

It is possible that both perceptions arise from a single underlying velocity field computation 
that combines velocity constraints along contours in the image, which is applied to image 
descriptions derived from initial V 2 G operators of different size. (Evidence for the existence of 
different size filters in early human vision can be found, for example, in Wilson and Bergen 
(1979).) Figures 16(b) and 16(c) show the zero-crossing contours derived from the convolution of 
the image of Figure 16(a) with V 2 G operators whose central positive diameter, w, is 6 and 16 
picture elements, respectively. With w = 6, the zero-crossing contours surround individual bars, 
while in the case where w = 16, the zero-crossings merge from one column to the next, forming 
continuous contours across the image. Figures 16(d) and 16(c) illustrate the initial perpendicular 
components of velocity, obtained from the temporal derivative and gradient of V 2 G*/. Figures 
16(f) and 16(g) illustrate the resulting smoothest velocity field in each ease. For the smaller 
operator, the velocity directions are roughly vertical, while for the larger operator, the smoothest 
velocity field indicates horizontal movement of the contours. Horizontal motion was perceived 
under the conditions of peripheral viewing or blurring, for which a coarser initial filtering of the 
image is expected (Wilson and Bergen 1979). 

5. Other Examples of Ideal Smooth Cur\’es 

Other examples of ideal smooth curves in motion, which illustrate the consistency between 
the velocity field of least variation, and human motion perception, shown in (Hildreth 1984a, b), 
include the kinetic depth effect stimulus, used by Wallach and O’Connell (1953) to demonstrate the 







a. b. 

Figure 17. A natural image sequence, (a) The original image, containing 576 x 576 picture elements, (b) 
The resulting zero-crossing contours.(The images in (a) and (b) are mirror reversed.) 

ability of the human visual system to derive three-dimensional structure from motion, the rotating 
deformed circle studied in (Walluch el al. 1956), and a rotating helix, which yields a perception 
of motion that is similar to the barberpole illusion (the perceived downward motion of the stripes 
of a barberpole). 


5.4.2 Natural Image Sequences 

Ihe final set of examples includes contours that have been extracted from natural images. 
For the first two examples, the motion was constructed artificially. That is, a single natural image 
was translated or rotated to produce a second image. In this case, the known velocity field allows 
a rigorous evaluation of the correctness of the resulting solution. The final example is a sequence 
of aerial photographs. Here, we rely on a qualitative assessment of the results. 


/. Pure translation of a curve from a natural image 

The image of Figure 17(a) was translated to the right by three picture elements, to yield 
a second image. The two images were convolved with a V 2 (7 operator, whose central positive 
diameter was 12 picture elements, and the zero-crossings, shown in Figure 17(b), were obtained. 
1Tic time derivative, £(V 2 <7*7) was computed by subtracting the two filtered images. At the 
location of a zero-crossing, u- 1 - was computed from the gradient of the filtered image, and u- 1 was 
computed using Equation (13). 

The analysis of the velocity field is presented for the zero-crossing contours extracted from 
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Figure 18. Translating contours, (a) Zero-crossing contours, shown with the true velocity field (b) The 
initial perpendicular components of velocity, (c) The computed velocity field. 
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the area outlined by the rectangle in Figure 17(a), The contours, together with die Lruc velocity 
field, arc shown in Figure 18(a). The initial perpendicular velocities arc shown in Figure 18(b). In 
__ Figure 19, the perpendicular vectors tAu- 1 along the contours arc displayed as points in velocity 

(fl 

space. If the initial measurements were perfect, the points in velocity space would lie along the 
solid circle shown. The deviation of the position of the points from this circle provides a visual 
demonstration of the large error in the initial motion measurements. The resulting velocity field 
solution, with /3 — 0.001, is shown in Figure 18(c). The average error in the direction of velocity 
is 2.2°, and average error in magnitude is 2.8% of the true magnitude of velocity. The directions 
of all resulting velocity vectors arc within 5° of the true direction of velocity. It is significant that 
a velocity field solution can be obtained, for which the error in velocity is so small, given the 
large error in the perpendicular components of velocity that formed the input to the algorithm. 
This error in the input underscores the need to design an algorithm that only requires the velocity 
field to satisfy the image constraints approximately. 

2. Rotation of a cur\e from a natural image 

For the next example, the original image was rotated rigidly to obtain a second image. 
Figure 20(a) shows the zero-crossing contour that is derived from the rectangle of Figure 17(a), 
and its true velocity field. Hie original V 2 C operator had a central positive diameter of 16 picture 
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Figure 20. Rotating contour, (a) A zero-crossing contour, shown with its true velocity field, (b) The initial 
perpendicular components of velocity, (c) The computed velocity field, (d) The true and computed velocity 
fields superimposed. 
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elements. 'ITic image was rotated .'5° in the clockwise direction, about a point in the center of the 
palm of llie hand. (The length of die displacement vectors in Figure 20 is slightly exaggerated.) 
Again, die initial perpendicular components of velocity were computed from the time derivative 
between the two filtered images, and die gradient along the /cro-crossing contours of the first 
filtered image: these components arc illustrated in Figure 20(b). Ihc computed smoothest velocity 
field, with 0 = 0.001, is shown in Figure 20(c). The true and computed velocity fields arc shown 
superimposed in Figure 20(d). The directions of velocity for the true and computed motions agree 
quite well around the fingers, thumb, and rightmost boundary of the hand. Around the palm of 
the hand, and the contour running from die thumb to the forefinger, there is considerable error 
in the direction of velocity. Error in the area of the palm of the hand may be due to a reduced 
contrast along the contour. The average error in the direction of velocity over the entire contour 
is 8.1 The magnitudes of velocity are generally smaller for the smoothest velocity field: average 
error in magnitude is 13.1% of the true magnitude of velocity. The solution varies with the choice 
of 0; a larger value of 0 leads to better agreement in the magnitudes of velocity between the true 
and smoothest velocity fields, but the error in direction of velocity increases. 

In general, the velocity field of least variation, for smooth curves in rotation, is expected to 
differ from the true velocity field. This is true for rotation in the image, as well as rotation in 
space, and was demonstrated for a number of ideal smooth curves in the previous section. The 
error in the velocity field of Figure 20 is due in part to error in the initial motion measurements, 
and in part to the fact that even with ideal input, we do not expect to obtain the correct velocity 
field in this case. 

3. A natural motion sequence 

Figures 21(a) and 21(b) show two aerial photographs, taken in sequence from an airplane. 
The images contain 256 x 256 picture elements. The two images were convolved with a V l G 
operator whose central diameter was 12 picture elements. The resulting zero-crossings, for the 
two images, are shown in Figures 22(a) and 22(b). The contours cover a large extent of the 
image in this case; this is a consequence of the small size of the image, and large displacements, 
which required a large initial operator. If the image were analyzed with higher spatial 
and temporal resolution, the two-dimensional velocity field could also be obtained with higher 
resolution. In Figure 23. the two zero-crossing descriptions arc shown superimposed, in order 
to illustrate the relative displacements between the two images. The zero-crossings from the first 
imafee are shown in white, those from the second are black, and the background is grey. Points 
at which zero-crossings occurred in both images are shown in white. Qualitatively, it can be seen 
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Figure 21. A natural image sequence, (a) and (b) Two natural images, containing 256 x 256 picture 
elements, taken in sequence from an airplane. 





Figure 21 Initial zero-crossing descriptions, derived from the images in Figure 21. 


that displacement increases in magnitude as we move from the top to the bottom of the image. In 
addition, the displacements vary in direction over the image, having a larger horizontal component 
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toward the top of the image, and a larger vertical component toward the bottom. 


As before, the initial perpendicular components of velocity were computed along the zero¬ 
crossing contours in the first filtered image by comparing the local gradient at the zero-crossings 
with the time derivative, obtained by subtracting the two filtered images. The algorithm was then 
run over each contour separately, to obtain the velocity field along the zero-crossing contours of 
the first image. In Figure 24, the zero-crossings are shown in black, with the resulting velocity 
vectors in white. The vectors shown represent a sampling of the velocity field along the contours. 
Very small contours, and some of the contours that occurred at the image boundary, were not 
included in this analysis. The length of the vectors in this display is equal to their displacement 
in the image. Qualitatively, the magnitude of the displacement vectors increases from the top to 
the bottom of the image, and their directions also vary in a way that reflects the displacement of 
the zero-crossing contours in Figure 23. 
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In Figure 1, shown in Section 2, the first image of this sequence, reduced in contrast, was 
displayed with the velocity field superimposed with black lines. The velocity vectors are evenly 
spaced on the image, and were computed by taking the average of the velocity vectors within a 
neighborhood of size 48 x 48 picture elements, centered at each point. Ihc length of the velocity 
vectors was doubled in this display, in order to emphasize the variation in direction and magnitide 
of velocity over the image. The resulting velocity field agrees qualitatively with the displacements 
of the contours shown in Figure 23. It also agrees well with the perceived movement between the 
images, when viewed in rapid succession. 

6. CONCLUSIONS 

In this paper, we have examined the computation of the two-dimensional velocity field 
that results from the changing projection of three-dimensional surfaces in motion. A theory for 
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this computation was developed, with three main components. First, initial measurements of 
motion in the image take place at the location of significant intensity changes, which give rise 
to zero-crossings in the output of the convolution of the image with a V'-Y; operator. The initial 
motion measurements provide the component of velocity in the direction perpendicular to the 
local orientation of the zero-crossing contours. Second, these initial measurements tire integrated 
along contours to compute the two-dimensional velocity field. Third, an additional constraint of 
smoothness of the velocity field, based on the physical constraint of smoothness of surfaces, allows 
the computation of a unique velocity field. 

There are three aspects of this work that are of fundamental importance. First, the additional 
constraint for the velocity field computation dial was formulated here allows the analysis of general 
motion, while providing a velocity field solution that is unique, and physically plausible. Second, 
this formulation of the computation leads naturally to algorithms that are biologically feasible, in 
that they require only simple, local, parallel operations. Finally, the computation appears to yield 
solutions that are consistent with human motion perception, and therefore may contribute toward 
our understanding of the measurement of motion in the human visual system. 
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AITINDIX A. I.YA1.HATING 011II R SM00TIINKSS Ml ASliRFS 

In this appendix, we first present an example for which the measures of variation given 
by /1‘^l^.s (variation in direction) and (variation in magnitude) do not yield unique 

solutions. We then present an example for which a velocity field that minimizes docs not 

exist. In both cases, the formulation of velocity constraints in velocity space is used to analyze the 
behavior of the smoothness measures. 

1. A Line Segment in Motion 

Consider a simple line segment that moves rigidly in the image plane, as shown in Figure 
25(a). l.ct us assume that the perpendicular component of velocity along the line, t>-L( fl ), and the 
velocities of the two endpoints, V(« 0 ) and V(s n ), arc known. The velocity of a particular point, 
V(.s t ) on the line is constrained to project to a line L, in velocity space, whose orientation is equal 
to the orientation of the line in the image, and whose perpendicular distance from the origin in 
velocity space is given by t>T(,s,). This constraint is illustrated in Figure 25(b). A set of evenly 
spaced points along the line yields a set of constraint lines in velocity space, which are evenly 
spaced from the origin, as shown in Figure 25(c). In addition, the velocity vectors associated with 
the endpoints are shown in velocity space. 

The true velocity field for the moving line segment corresponds to the straight line / in Figure 
26(a), connecting the two known velocity vectors. The velocity of a point on the line segment, 
V(a,), is given by the intersection of / and U. Any valid velocity field must correspond to a path 
in velocity space that begins at the endpoint of V(« 0 ), intersects each constraint line in the order 
that they arise along the image curve, and ends at the endpoint of V(« n ). Three examples of valid 
velocity fields correspond to the curves shown in Figure 26(b). Computation of the velocity field 
in this case involves finding the valid velocity field that exhibits the least variation, given one of 
the measures of variation. 

Suppose that we want to find die velocity field that is consistent with the constraints shown 
in Figure 25(c), and minimizes the total change in direction of velocity, given by: / |§f | da. For 
the case of the true velocity field, ^ increases monotonically; its total change is given by tine angle 
«, shown in Figure 26(a): hence / |§£j ds — n. n is the smallest total variation in direction that a 
velocity field consistent with these constraints may obtain. There are other velocity fields, however, 
that distribute the total change « over the curve monotonically and smoothly, in a way that is 
consistent with i>-J-(,s). Other examples correspond to the curved paths in Figure 26(b); for these 
velocity fields, it is also the case that /1 \ds — rv. Hence, in the simple case of a line segment 













Figure 25. A line segment in motion, (a) A line segment moving rigidly in the image, with known velocity 
vectors at its endpoints, (b) The constraint imposed by a single measurement of t>- L (s) (shown in the image 
and in velocity space), (c) The constraints imposed by multiple measurements of e- L (s) along the line 
segment 
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Figure 26. Velocity fields for a line segment in motion, (a) The true velocity field in the image projects 
to the line / in velocity space, (b) The curved paths correspond to other valid velocity field solutions. 

moving rigidly in the image plane, minimizing the total change in direction of velocity alone docs 
not yield a unique velocity field solution. 

Now consider minimizing the variation in magnitude of velocity, given by the measure: 
/ Let and m 2 denote the magnitudes of the velocity vectors associated with the 

endpoints. It then follows that / > |mi — m 2 \. Any velocity field that is consistent 

with the image constraints, and changes monotonically in magnitude along the line will yield 
= |mi - m 2 |. This solution is again not unique. 


2. A Rotating Semi-circle 

As a second example, consider a semi-circle rotating about its center, as shown in Figure 
27(a). Assume that the velocity vectors V(s 0 ) and V(«„) at the endpoints are known. Along the 
contour, v-L( s ) = o, so these measurements give rise to lines of constraint in velocity space that 
intersect at the origin, as shown in Figure 27(b). The known velocity vectors at the endpoints are 
also shown in Figure 27(b). Consider the measure of variation given by: f\^j\ds. As before, a 
valid velocity field must correspond to a path in velocity space that begins at the endpoint of V(« 0 ), 
intersects the lines of constraint in the order that they arise from the image curve, and ends at the 
endpoint of V(a„). The value of the measure is equal to the length of the path in velocity 

space. Thus, a velocity field that minimizes this measure must give rise to a path of minimum 
length in velocity space. For this particular example, however, a minimum length path satisfying 
the constraints docs not exist. In Figure 27(c), we show three paths in velocity space, labeled 
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P u P-i and P 3 , that correspond to valid velocity fields. As we move from Pi to P 3 , the length 
of the path decreases. Tims, the value of the measure J\^j\ds decreases, for the corresponding 
velocity fields. We can continue indefinitely, however, to find paths of shorter length that satisfy 
the constraints. The paths will approach the horizontal line, but can never reach it exactly, because 
the horizontal line docs not correspond to a valid velocity field. Thus, a velocity field of least 


variation, given this particular measure of variation, does not exist. 



APPENDIX 15. TUI CONJUGATE GRADIENT ALGORITHM 


In Section 5, it was shown that the velocity field computation, in which the constraints 
imposed by the measurements of o 1 (,s) arc only approximately satisfied, can be formulated as an 
unconstrained optimization problem, and the conjugate gradient algorithm can be used to obtain 
a solution. The computation finds the x and y components of velocity for each of the » points on 
a contour. Each possible solution can be considered as a point in a 2u-dimensional vector space, 
with each dimension representing the x or y component of velocity at one point on the contour. 
The 2n-dimcnsional vector space can be embedded in a (2 n + l)-dimcnsional vector space, w'hcre 
the final dimension corresponds to the total variation of velocity along the contour. By associating 
a measure of variation with each possible velocity field, we can construct a hypersurface , called 
the objective surface, in the (2n + l)-dimensional space. The goal of the velocity field computation 
is to find the "lowest" point on this surface; this point corresponds to the velocity field of least 
variation. 

The conjugate gradient algorithm is an iterative descent algorithm; a sequence of approximations 
V (1) ,V< 2 ),... to the exact solution V is computed, given an initial approximation V ( °l, and this 
sequence has the property that each new approximation decreases the value of the objective 
function. The basic steps for a descent method are as follows: 

(1) Start at an initial point V { °1, which is usually V(°> = 0. 

(2) According to a fixed rule, determine a direction of movement along the objective 

surface that will reduce the value of the objective function. 

(3) Move in this direction to a (relative) minimum of the objective function. 

(4) If the final solution has not yet been reached, return to step 2. 

One of die most common descent methods is that of steepest descent, in which the direction of 
movement in step 2 above is given by the negative gradient of the objective function, evaluated at 
the present point. The steepest descent method has the disadvantage of relatively slow convergence. 
A much improved rate of convergence can be obtained with the method of conjugate gradients. 

The velocity field of least variation is given by the solution of a system of linear equations. 
Let diis system be denoted in matrix form by: 

Qx — b 

where Q is an n x n symmetric positive definite matrix. The unique solution to diis equation is 
equivalent to die unique solution of the quadratic problem: 








minimize -x r Qx-b T x 
subject to x G n C E n . 


l'J n denotes the n-dimensional Huclidcan space. Conjugate direction methods were invented largely 
for die analysis of quadratic problems. We present a brief development of these methods, and the 
conjugate gradient algorithm. This development is taken from I.uenbergcr (1973), where further 
details may be found. 

The basic idea behind conjugate direction mcdiods is to optimize die direction in which 
each new step is taken in the algorithm. To show how this optimization is achieved, we begin by 
defining the notion of Q-orthogonality: 


Definition. Given a symmetric matrix Q, two vectors di and d 2 are said to be Q-orlhogonal, 
or conjugate, with respect to Q if d7Qd 2 = 0. 

The following can then be shown: 


Proposition. If Q is positive definite and the set of nonzero vectors d 0 , di,...,d* are 
Q-orlhogonal, then these vectors are linearly independent. 

Given the positive definite matrix Q, and a set of non-zero Q-orthogonal vectors, d 0 , dj, ...d„_i, 
their linear independence allows us to express the solution x' in terms of these vectors: 


x' — C*odo + ••+ a n-ld n -l 




for some set of scalars a,-. The orthogonality of die vectors implies that multiplying by Q and 
taking the scalar product with di causes all the terms, except the i th to vanish, yielding: 

d T Qx' d t T b 
Q ‘“dTQd,~d7Qdi- 


The solution is then given by the following: 


V 1 d C . b d- 
toAjQii 


The expansion for the solution x' can be considered as the result of an iterative process of n steps, 
where at the i th step, the term a,d, : is added. In this way, die following result can be proven 
(Lucnberger 1973): 
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Conjugat e Direction Theo rem: Let {d,}d be a set of non-zero Q-ortliogonal vectors, 
lor any x 0 E / i" the sequence {x t } generated according to 


x*+ i = X*: + fVfcdfc , k <0 


with 


and 


«k 


gjd* 




gk = Qxk - b 

converges to the unique solution. x\ of Qx = b after n steps. 

The conjugate gradient algorithm provides one method for generating the sequence of 
Q-orthogonal direction vectors; the successive direction vectors are selected as a conjugate version 
of the successive gradients obtained as the method proceeds. Hence, at each step, the next direction 
vector is determined, based on the current state of the objective function and its gradient vector. 
The algorithm, which is guaranteed to converge in at most n steps, is outlined below: 

The Conjugate Gradient Algorithm 

Starting at any point x 0 E E n define d 0 = -go = b - Qx 0 and 

Xk+i = x* + «fcd* 

6k<U 

* d£Qd* 

dk+l = -gk+l + Pk&k 


/9k 


gk+lQdk 


djQdk 

where g* = Qx* - b. 

For the velocity field computation, we showed in Section 5.2 that the system of 2 n linear 
equations to be solved is given by: 


3<l>o 


0, 




a\ Xi ’ d\y 


o, 


1 < i < n 


where $ 2 is given by: 

*1 = E [O'- - v„„) ! + (V„ - V„_,) ! ] + [(V., - V,. )* + (V„ - V,.)«] 


1=2 
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fT. 


i=i L 


t2 


4 \’v,uj:-v± 


■rom the partial derivatives with respect to die x components of velocity, we have the equations: 
(4 + ‘W(»±f)\’ Xi - 2V Ij+I - 2V*,_, 4 2f)a^\ y , = 2/3v±u± it 1 <i<n. 

’rom the partial derivatives with respect to the y components of velocity, we have: 


(4 4 2P(u±)*)\ v , - 2V„ +1 - 2V y> _, 4 2flu±u±\ x , = 2(h^u±, 1 <i<, 


The coefficients on the left hand side of the equations form the symmetric positive definite matrix 
Q, and the values on the right hand side form the column vector b. It should be noted that for 
the case of an extended straight line, Q is singular, so that a unique solution does not exist. The 
initial velocity field consists of the normal velocity vectors; hence: 


Xo = (»f ....^u^wfujj^u. 


J-ii-L „-Ln_L 


> V n U ^) 


The conjugate gradient algorithm is then applied straightforwardly. 
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